Polynomial-time quantum algorithm for the simulation of chemical dynamics 
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The computational cost of exact methods for quantum simulation using classical computers grows exponen- 
tially with system size. As a consequence, these techniques can only be applied to small systems. By contrast, 
we demonstrate that quantum computers could exactly simulate chemical reactions in polynomial time. Our 
algorithm uses the split-operator approach and explicitly simulates all electron-nuclear and inter-electronic in- 
teractions in quadratic time. Surprisingly, this treatment is not only more accurate than the Born-Oppenheimer 
approximation, but faster and more efficient as well, for all reactions with more than about four atoms. This 
is the case even though the entire electronic wavefunction is propagated on a grid with appropriately short 
timesteps. Although the preparation and measurement of arbitrary states on a quantum computer is inefficient, 
here we demonstrate how to prepare states of chemical interest efficiently. We also show how to efficiently 
obtain chemically relevant observables, such as state-to-state transition probabilities and thermal reaction rates. 
Quantum computers using these techniques could outperform current classical computers with one hundred 
qubits. 



Accurate simulations of quantum-mechanical processes 
have greatly expanded our understanding of the funda- 
mentals of chemical reaction dynamics. In particular, re- 
cent years have seen tremendous progress in methods de- 
velopment, which has enabled simulations of increasingly 
complex quantum systems. While it is strictly speaking 
true that exact quantum simulation requkes resources that 
scale exponentially with system size, several techniques 
are available that can treat realistic chemical problems, at 
a given accuracy, with only a polynomial cost. Certain 
fully quantum methods — such as multiconfigurational time- 
dependent Hartree (MCTDH) matching pursuit/split- 
operator Fourier transform (MP/SOFT) [2], or full multiple 
spawning (FMS) [,3J] — solve the nuclear Schrodinger equa- 
tion, including nonadiabatic effects, given analytic expres- 
sions for the potential energy surfaces and the couplings be- 
tween them. These techniques have permitted the simulation 
of large systems; as examples we can give MCTDH simula- 
tions of a penta-atomic chemical reactionM and of a spin- 
boson model with 80 degrees of freedom ||5|], or an MP/SOFT 
simulation of photoisomerization in rhodopsin using 25 de- 
grees of freedom [6]. Ab initio molecular dynamics tech- 
niques such as ab initio multiple spawning (AIMS) ||7|] avoid 
analytic expresions for potential energy surfaces and instead 
solve electronic Schrodinger equation at every timestep. This 
allows one to gain insight into dynamical problems such as 
isomerizations through conical intersections |8]. 

However, there are also chemical processes which are best 
treated by completely avoiding the Born-Oppenheimer ap- 
proximation. As examples we can cite strong-field electronic 
dynamics in atoms and multi-electron ionization or 
atomic and molecular fragmentation caused by collisions with 
energetic electrons or photons il liil2ll . Systems that resist the 
application of the Born-Oppenheimer approximation require 
very general techniques, and the consequent unfavorable scal- 
ing has restricted such simulations to systems with a few parti- 
cles. Here, however, we show that the Born-Oppenheimer ap- 



proximation would not necessarily simplify simulations per- 
formed on quantum computers. Indeed, except for the small- 
est systems, an explicit treatment of all the particles would be 
both more accurate and more efficient, even for nearly adia- 
batic chemical reactions. 

Feynman's idea of using a quantum machine to mimic the 
quantum Hamiltonian of a system of interest was one of the 
founding ideas of the field of quantum computation 111 311 . 
Lloyd [14] subsequently showed that quantum computers 
could be used to simulate systems which can be formulated 
in terms of local interactions, using resources that scale only 
polynomially with system size. Zalka and Wiesner fisl [3 
developed a quantum simulation algorithm for particles in real 
space and Lidar and Wang IitIi applied it to the calculation of 
the thermal rate constant of chemical reactions. Smirnov et 
at. ifisl ] proposed an analog quantum simulator for chemical 
reactions using quantum dots. We have previously shown lll9ll 
that quantum computers could be used to simulate the static 
properties of molecules, and in this work we present a general 
scheme for using quantum computers for the study of dynam- 
ical chemical properties. 

To simulate a quantum system we must prepare its initial 
quantum state, propagate it in time, and finally extract data 
of chemical relevance, such as rate constants. For an efficient 
quantum simulation, all these tasks must be carried out using 
resources which increase polynomially with increasing sys- 
tem size. We present a quantum algorithm that meets these re- 
quirements. We also show that for all chemical reactions with 
more than about four atoms, it is more efficient for a quan- 
tum computer to simulate the complete nuclear and electronic 
time-evolution rather than to use the Born-Oppenheimer ap- 
proximation. 

The polynomial scaling of these methods means they would 
enable the study of systems which are in principle out of reach 
for any classical computer. However, large quantum comput- 
ers are far in the future, and so determining the requirements 
of interesting calculations in absolute terms is, perhaps, of 
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more interest than their scaling alone. We show that a quan- 
tum computer using these techniques could outperform cur- 
rent classical computers using one hundred qubits, within the 
design limits of a proposed 300-qubit quantum computer i20ll . 
While we focus on chemical applications, these techniques are 
generally applicable to many physical systems, from strong- 
field, multielectron ionization to quantum liquids and con- 
densed matter systems. 

This article is organized as follows. We first review Zalka 
and Wiesner's algorithm and show how the difficulty of com- 
puting a wavefunction's time evolution depends only on the 
complexity of evaluating the interaction potential. We then 
consider three approaches to the calculation of the interaction 
potential, including a fully non-adiabatic treatment of chem- 
ical reactions. We consider the problem of state preparation 
for all of the schemes, and finally address the problem of 
measurement. We present three readout schemes for reaction 
dynamics — ^reaction probabilities, thermal rate constants, and 
state-to-state probabilities — which would allow for efficient 
evaluation of many parameters accessible to experiment. 



QUANTUM DYNAMICS 

The problem of simulating quantum dynamics is that of de- 
termining the properties of the wavefunction \ip{t)) of a sys- 
tem at time t, given the initial wavefunction \ip{0)) and the 
Hamiltonian H of the system. If the final state can be prepared 
by propagating the initial state, any observable of interest may 
be computed. 

We employ an improved version of the real-space quan- 
tum simulation technique developed by Zalka and Wiesner in 
which a discrete variable representation of the wavefunction 
is used 115, 3. In the one-dimensional case, the domain of 
the wavefunction is divided into a discrete position basis of 
N = 2" equidistant points. The wavefunction is represented 
as: 
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The spatial wavefunction is stored in the Hilbert space of the 
qubits, and so the spatial resolution grows exponentially with 
the number of qubits. For a system with d dimensions, d reg- 
isters of n qubits each are used, |x) = • • • la;^), represent- 
ing a grid of 2*^" points. The states of multiple particles can 
be stored by adding position registers for each particle. There- 
fore, only a polynomially large number of qubits is required 
to store the system wavefunction. 

For simplicity we assume a time-independent Hamiltonian 
whose potential depends only on position, H = T + V 
where T = p^/2m and V = V{x) are the kinetic and 
potential en ergy op erators, respectively. The split operator 
method lllSi l2ll 12211 computes the time evolution by separat- 
ing the kinetic T and potential V energy contributions to the 
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FIG. 1: The quantum simulation algorithm. The potential and ki- 
netic energy unitaries are applied to a quantum state in turn, with the 
transformation between position and momentum representations be- 
ing performed with the efficient quantum Fourier transform (QFT). 
The ancilla register is required for phase kickback and remains un- 
changed throughout the simulation, while the boxed time step is re- 
peated t/St times. The proposed algorithm, unlike that of Zalka 1 1 5] , 
does not require that functions be uncomputed and is therefore twice 
as fast. 



propagator U{t) = e Given a sufficiently small time 

step St, we can write to first order 
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The operators e~*^''* and e~'^'^^* are diagonal in the posi- 
tion and momentum representations, respectively. A quantum 
computer can efficiently transform between the two represen- 
tations using the quantum Fourier transform (QFT) |23]: 
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The procedure is iterated as many times as necessary to obtain 
the system wavefunction \ip{t)) after an arbitrary time t to a 
desired accuracy. 



The application of diagonal unitaries is straightforward 
on a quantum computer Suppose that we have a gate se- 
quence which acts on an arbitrary position eigenstate as |x) 
^-zV{x)St Since lip) is a superposition of position eigen- 
states, when this gate sequence is applied to one obtains 
^-iVSt ^ single application. 

We depart from Zalka and Wiesner's method in the imple- 
mentation of this gate sequence. We are free to take the lowest 
value of the potential in the domain as 0, and use such units 
that the maximum value of the potential is Knax — 2™ — 1, 
with m an integer. With this choice of units V takes inte- 
ger values, and we choose m large enough that V is resolved 
with sufficient precision. The integer m is therefore the num- 
ber of qubits required to represent the desired range of po- 
tential values with the desired precision. The gate sequence 
V which computes the potential V acts so that V|x, y) = 
\x,y ® V{x.)), where y is an m-bit integer labeling a basis 
state of the ancilla register and denotes addition modulo 
2™. 

We apply the diagonal unitary by phase kickback. The com- 
puter is initialized in the state \tp) |l)„j, where \ l)„^ in the 
ancilla register represents the state |0. . .001) in ni qubits. 
Applying the inverse QFT to the ancilla register, followed by 
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V, produces 
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where M — 2™ and we choose 5t — The equal- 

ity obtains since the ancilla state is an eigenstate of addi- 
tion (with eigenvalue e^'^^'^i/^'^ corresponding to the addi- 
tion of q) 1 ,24] , We see that applying V results in the req- 
uisite diagonal unitary action on the wavefunction register. 
The states of register and ancilla are separable before and af- 
ter each potential evaluation. We can also define a quantum 
gate sequence T which computes the kinetic energy /2m: 
^ |P> v) — IPi y ffi ^(p))- This gate is diagonal in the mo- 
mentum basis, and has efficiently computable entries on the 
diagonal (namely p^). Thus, we use the quantum Fourier 
transform to conjugate into the momentum basis and T is im- 
plemented by phase kickback in exactly the same way as V. 
The quantum circuit for this algorithm is shown in Fig.[T] 

This simulation algorithm is numerically exact in the sense 
that all introduced approximations are controlled, so that the 
error in the calcuation can be arbitrarily reduced with an ad- 
ditional polynomial cost. The only approximations employed 
are the discretization of time, space, and the potential V^(x). 
The error due to discretization can be made exponentially 
small by adding more qubits. The error due to time discretiza- 
tion can be systematically improved by use of higher-order 
Trotter schemes 1I25I1 . The computational cost of the algo- 
rithm per iteration is the evaluation of VipC), T(p) and two 
QFT's. While the QFT's and the quadratic form in the kinetic 
energy {p^ in the simplest case) can be computed in polyno- 
mial time ['23', '2^, the evaluation of the potential energy F(x) 
may not be efficient in general. For example, a random poten- 
tial stored in an exponentially large database requires an ex- 
ponentially large number of queries to evaluate. However, any 
classical algorithm running in 0{f{n)) time can be adapted to 
a reversible quantum algorithm also running in 0{f{n)) time 
||27|] . Therefore, the potential energy ^(x) will be efficiently 
calculable on a quantum computer if it is efficiently calcula- 
ble on a classical computer. Fortunately, this is true for all 
chemically relevant cases. 



CHEMICAL DYNAMICS 

Every isolated system of chemical interest has the same 
Hamiltonian, which in atomic units is 



H 



mi 



where the sums are over the nuclei and electrons, pi is the mo- 
mentum of the i*'' particle. Mi its mass, qi its charge and 
is the distance between particles i and j. Both the potential 
and kinetic terms can be efficiently evaluated since the arith- 
metical operations can be performed in 0{m?) time ll26ll . and 



for a system of B particles, there are only 0{B'^) terms in the 
sum. |43!l 

The fact that the Coulomb potential can be evaluated in 
0{B^rin?) time implies that chemical dynamics could be sim- 
ulated on a quantum computer in 0{B^m^) time, an expo- 
nential advantage over known classical techniques for exact 
quantum simulation. Here B is the number of particles and m 
is the binary precision the potential in the region of interest. 
We want to emphasize that a quantum simulation would be 
substantially different from what is usually done on classical 
computers. Most significantly, we are proposing to explicitly 
track all the nuclei and electrons on a Cartesian grid which is 
sufficiently fine and with time steps sufficiently short to cap- 
ture the true electronic dynamics. We will show that this is 
not only more accurate, but also requires fewer quantum re- 
sources. 

The Supplementary Information contains a detailed com- 
putation of the numbers of gates and qubits required for the 
quantum simulation of the Coulomb potential. The number 
of elementary gates required to evaluate this potential in three 



dimensions is (-j-m 



51, 



per pair of particles (Fig. 



We chose a method which minimizes the number of ancilla 
qubits and so is suited for small numbers of qubits. Note that 
this scaling is not asymptotically optimal (the asymptotic re- 
quirement would be 0{m^)), so further improvement could 
be achieved for computations with high precision (large m) 
if suitable arithmetical algorithms were implemented. Storing 
the wavefunction of a system with d degrees of freedom re- 
quires nd qubits, so a system of B particles, with d = 3B — 6 
degrees of freedom, requires n{3B — 6) qubits. To this one 
must add the qubits needed for the ancilla registers, only four 
of which are required for the Coulomb potential, meaning that 
simulating these potentials requires n{3B — 6) + 4m qubits 
(Fig.EJ. 



On a small quantum computer, the computational cost of 
simulating the interactions between many particles may prove 
prohibitive. One could try to simplify matters and focus only 
on the nuclear motion by employing the Born-Oppenheimer 
approximation. Here, the solution of the electronic struc- 
ture problem provides a potential energy surface on which 
the nuclei move. However, we show that quantum comput- 
ers would benefit from the Born-Oppenheimer approximation 
only rarely: for many chemical reactions, simulating the full 
dynamics of both electrons and nuclei will not only yield more 
accurate results, but will also, remarkably, be faster. This is in 
sharp contrast to the study of chemical dynamics on classical 
computers, where there is frequent need to simplify calcula- 
tions using the Born-Oppenheimer approximation. 

It is difficult to estimate the precise computational cost 
of using the Born-Oppenheimer approximation, since differ- 
ent potential energy surfaces have different functional forms. 
Nevertheless, for any general fitting technique, the complex- 
ity of the interpolating function grows exponentially with in- 
creasing dimension of the system, since exponentially many 
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FIG. 2: Resource requirements for a quantum simulation of B par- 
ticles interacting through a pairwise potential. The chemical sym- 
bols correspond to the simulation of the full Coulomb dynamics of 
the corresponding atom (nucleus and electrons). The vertical dashed 
line represents the approximate current limit of numerically exact 
quantum simulation on classical computers on a grid fioll . (A) Total 
qubits required. We require n qubits for each degree of freedom and 
m qubits for each ancilla, four of which are needed for the Coulomb 
potential. Hence, a total of n{3B — 6) + 4m qubits are needed (see 
Supplementary Information for details). The horizontal dotted line 
represents a 300-qubit quantum computer, which is believed to be 
feasible with near-future technology ||2(?]. We assume a grid of 2'"' 
points, which corresponds to n = 10 and would suffice for the sim- 
ulation of many chemical reactions or the strong-field ionization of 
atoms |9, 28]. (B) Total elementary gates required. The 300-qubit 
computer is expected to achieve one billion elementary quantum op- 
erations. The dotted line represents the largest possible simulation of 
1000 time steps, assuming ten bits of numerical accuracy (m — 10). 
Computing the Coulomb potential requires {■^m'^ + ^-ttl^) gates 
per pair of particles (see Supplementary Information for details). 



data points must be used in the fit if one is to maintain uni- 
form accuracy over the entire domain. We can provide an es- 
timate of the computational cost for a potential energy surface 
which is an interpolating polynomial of degree K along each 
dimension (see Supplementary Information). In that case, the 
total cost of the adiabatic simulation is K'^^^^ (fm-"^ + 1"^-^) 
per nuclear time step (which is usually about a thousand times 
longer than an electronic time step). Numerical experiments 
with the BKMP2 surface for H3 indicate that K must be 
chosen to equal at least 15 if one aspires to 0.1% accuracy in 



the potential, and more for chemical accuracy. With K = 15, 
the exponential growth implies that even for heavy elements 
(Z sa 100), the fully-dimensional diabatic treatment is faster 
for all reactions involving more than four atoms, and even for 
many smaller reactions, as shown in Fig. [3] 

It is perhaps beneficial to briefly discuss the intuitive rea- 
sons why the use of pre-computed potential energy surfaces 
is not as useful on quantum computers as it is on classical 
machines. Classically, an exponential amount of memory is 
required in general to store the wavefunction. However, the 
ratio of computing power to memory in most classical com- 
puters is large and the basic floating point operations are hard- 
wired. Since the storage capacity is often the limiting factor, 
if a wavefunction can be stored in memory, its motion on a 
surface can probably be computed. Quantum computers, on 
the other hand, require only linearly many qubits to store the 
wavefunction in a superposition state. However, using a pre- 
computed potential requires either the evaluation of a com- 
plicated function or a look-up in a potentially large table. The 
potential energies must be computed on the fly in order to take 
advantage of quantum parallelism and it is therefore imper- 
ative to keep the interaction potential as simple as possible. 
This is achieved by treating all the particles explicitly, inter- 
acting via the Coulomb interaction. 

An alternative way to compute a potential energy surface 
would be to embed an on-the-fly calculation of electronic 
structure in the quantum algorithm and thus avoid a classi- 
cally precomputed fit. This can be done efficiently since elec- 
tronic structure calculations can be performed in polynomial 
time on quantum computers lll9ll . Hence, the quantum cir- 
cuit V would be replaced by a black box containing the ef- 
ficient quantum version of the full configuration interaction 
(FCI) procedure [19]. Because the quantum simulation algo- 
rithm exploits quantum effects, a single evaluation of elec- 
tronic structure is sufficient for each time step: all the nuclear 
configurations are evaluated in superposition. However, the 
electronic structure circuit for the proposed algorithm would 
require the atomic positions as input. This would require a 
modification of the original algorithm so that the Coulomb 
and exchange integrals are computed using a quantum circuit 
rather than classically. While this approach, unlike the Born- 
Oppenheimer approximation, is asymptotically appealing, the 
large overhead required to compute the exchange integrals 
quantum mechanically makes it uninteresting for near-future 
implementation. 

Steane has recently proposed a design for a 300-qubit, 
error-corrected, trapped-ion quantum computer that could 
perform around 10^ quantum operations using methods for 
quantum gates that have already been experimentally imple- 
mented 1 20]. On a three-dimensional grid of 2"^" points, such 
a computer could store the wavefunction of a ten-particle sys- 
tem (Fig. |2]i. By comparison, classical computers implement- 
ing a comparable grid based algorithm are limited to com- 
puting the full quantum evolution of a three-particle system, 
such as a helium atom [9, 10]. Even a relatively modest quan- 
tum computer with 100 qubits could simulate the electron dy- 
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FIG. 3: Estimated number of elementary quantum operations (gates) 
required for the simulation of chemical reactions. Standard Born- 
Oppenheimer potential-energy-surface calculations require time re- 
sources exponential in the size of the system (full line), while an fully 
nonadiabatic, nuclear and electronic calculation would require only 
polynomial time (dotted). The resulting cutoff indicates that for all 
reactions with more than four atoms (dashed) the Born-Oppenheimer 
approximation is always less efficient on a quantum computer than 
a diabatic treatment. The complexity of the diabatic computation 
depends only on the atomic number Z, while the potential energy 
surfaces are modeled as polynomials of degree K along each axis. A 
value of > 15 is required to obtain 0.1% agreement with surfaces 
such as BKMP2 |29]. The position of the cutoff does not signifi- 
cantly depend on the accuracy of the evaluated potential (m). To 
obtain the gate counts, we assume 20 bits of accuracy (m = 20), 
enough for chemical precision. The gate counts reflect the fact that 
an appropriate nuclear time step is about 1000 times longer than an 
electronic time step. 



namics or ionization of the lithium atom, a task beyond the 
reach of classical computers using grid based methods |10]. 
The simplest chemical reaction, H + H2 ^ H2 + H, is a six- 
particle system, and could therefore be simulated by Steane's 
computer in a fully-dimensional diabatic regime. While other 
classical methods may be able to reach somewhat larger ex- 
amples, the exponential scaling of all known classical exact 
methods means that the examples given here are close to the 
crossover point between classical and quantum computing for 
chemical dynamics. There remain two questions: how to pre- 
pare the initial state of the quantum computer, and how to 
extract useful information out of the final state. 



nuclear and electronic wavefunctions, each in its own quan- 
tum register 

Nuclear motions can be expressed in normal mode coordi- 
nates if the displacements from equilibrium are small, which 
is the case in molecules at chemically relevant temperatures. 
The nuclear wavefunction is then, along each coordinate, a 
superposition of harmonic oscillator eigenstates, which are 
themselves products of Gaussians and Hermite polynomials. 
It is known that superpositions corresponding to efficiently in- 
tegrable functions can be prepared on a quantum computer in 
polynomial time ll5l bill . Therefore, Gaussian wavepackets 
and Hermite polynomials are efficiently preparable, meaning 
that we can prepare good approximations to molecular vibra- 
tional states and Gaussian wave packets. 

Gaussian wavepackets can also be used to prepare the 
electronic wavefunction. Indeed, it is customary in elec- 
tronic structure theory to expand molecular orbitals in terms 
of atomic orbitals, which themselves are superpositions of 
Gaussian-type orbitals. The orbital occupation numbers can 
be obtained from electronic structure calculations, including 
our earlier quantum algorithm lll9ll . Consequently, the occu- 
pied orbitals, which are superpositions of Gaussians, can all 
be prepared efficiently. 

One final consideration is the exchange symmetry of multi- 
electron wavefunctions. Abrams and Lloyd proposed a 
method for preparing antisymmetric wavefunctions Ii32,1 start- 
ing with a Hartree product of molecular orbitals. We propose 
to use this method for preparation of multi-electron wavefunc- 
tions, noting that it suffices to prepare the initial state with the 
coiTect exchange symmetry, since the exchange operator com- 
mutes with the Hamiltonian. 

Of course, other strategies for state preparation can be pur- 
sued, such as the phase estimation algorithm 113 311 . If we 
are able to prepare a state \S) that has a significant over- 
lap {S\E) with an eigenstate \E) (not necessarily the ground 
state), phase estimation followed by measurement will col- 
lapse the wavefunction to the desired eigenstate with proba- 
bility Alternatively, the ground state can be pre- 
pared by the adiabatic state preparation algorithm El. This 
is of particular significance to the simulation of full chemical 
dynamics, since the electronic ground state is usually a good 
approximation for the reactants. 



STATE PREPARATION 

The preparation of an arbitrary quantum state is expo- 
nentially hard in general Il30ll . Nevertheless, we show that 
most commonly used chemical wavefunctions can be pre- 
pared efficiently. Since the significant deviations from Born- 
Oppenheimer behavior occur during evolution and usually do 
not concern initial states, we will prepare the initial state us- 
ing the Born-Oppenheimer approximation. That is, the system 
wavefunction will be a product state I?/;) = IV'nuc) IV'cicc) of 



MEASUREMENT 

After preparing the initial state and simulating its time evo- 
lution using the methods described above, we must extract 
chemically relevant information from the final system wave- 
function. In general, quantum tomography is the most general 
approach to the estimation of an unknown quantum state or 
a quantum process 112711 by measuring the expectation values 
of a complete set of observables on an ensemble of identi- 
cal quantum systems. However, this full characterization of 
quantum systems always requires resources that grow expo- 
nentially with the size of the system. In order to avoid such 
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problems, alternative approaches for the direct estimation of 
certain properties of quantum dynamical systems have been 
recently developed 11341 13511 . Here we likewise show that the 
data of greatest chemical significance can be obtained directly 
with only a polynomial number of measurements. In particu- 
lar, we present algorithms for obtaining the reaction probabil- 
ity, the rate constant, and state-to-state transition probabilities. 

The reaction probability, given a certain initial wavefunc- 
tion of the reactants, is the likelihood of observing, after a 
sufficiently long time, the products of the chemical reaction. 
To find it, we divide the real-space domain of the wavefunc- 
tion into r disjoint regions corresponding to sub-domains of 
chemical interest. In chemistry these regions are typically a 
few simple polytopes. The simplest division is into only two 
regions, one for the reactants and one for the products, sepa- 
rated by the transition state dividing surface (TSDS). The re- 
action probability is the sum of the probabilities of finding 
the final wavepacket in the product region(s). It is straightfor- 
ward to construct a classical point location circuit for the func- 
tion -R(x) which, given a nuclear position vector x, identifies 
which region it is in by returning an integer label correspond- 
ing to that region. There is a corresponding reversible circuit 
that performs the transformation |x) \y) |x) \y ® -R(x)). 
We can apply this circuit to the final state \tp) = '^^ 1^)' 
to which we add an additional ancilla register with [log2 r~\ 
qubits. That is, applying this reversible circuit to j?/;) |0) yields 
X^x I-") Measuring the ancilla register will return 

i with probability Pi, which equals the probability of finding 
the wavepacket in the region i. We can obtain all the probabil- 
ities Pi by employing an ensemble measurement of the ancilla 
register. Since individual measurements are uncorrected, the 
error of the estimate of the probabilities decreases as l/VAf 
for M repetitions of the experiment. However, it is possible 
to achieve a convergence of which is the ultimate limit 
of quantum metrology [36], using techniques such as ampli- 
tude estimation 11371 13811 . Next, we use these disjoint regions 
to compute the rate constant. 

The rate constant k{T) at temperature T is a thermally 
weighted average of cumulative reaction probabilities ifjoll : 

1 f°° 



1 



2iThQ{T) 



Y,Pr{C,E)e-^/''-^AE, 



where E is the energy, Q{T) is the reactant partition function, 
and N{E) is the cumulative reaction probability, N{E) = 
-Pr(C) E)- The vector ^ is a specification of all the quan- 
tum numbers of the reactants and Pr{C, E) is the reaction 
probability given that the reactants are in the eigenstate speci- 
fied by C and E. The sum ranges over all possible ^, and with 
E from zero to a cutoff. Note that on a quantum computer 
the cutoff can be made exponentially large, or the energy step 
AE exponentially small, by adding more qubits. 

We can compute the rate constant on a quantum com- 
puter if we propagate in time not a pure wavefunction but 



the correct thermal mixed state. In that case, the expecta- 
tion value of the reaction probability would equal the rate 
constant, up to a known factor The required initial state 
is p(0) = J:^ ^ r {E, Tf 100 (C, E)) (00 (C, E)\ where 
T{E,T) ^ {exp{~E/kBT)AE/2TrhQ{T)Y^^ is the square 
root of the appropriate Boltzmann factor, C is a normaliza- 
tion constant, and j^o {Ci E)) is a real-space reactant eigen- 
function corresponding to quantum numbers ^ and energy 
E. If we propagate p{Q) for time t using the simulation 
algorithm, the system will be in the final state p{t) — 
ECE r (E, Tf \cbt iC, E)) (C, E)\, where \cbt (C, E)) 
now denotes the time-evolved version of state \(j)Q{C,E)) 
(note that, except in exceptional cases, |0t {C, E)) is not an 
eigenfunction of either reactants or products). If we have a 
quantum register in this mixed state, we can add an ancilla 
qubit and use the technique of dividing the domain into reac- 
tant and product regions as described above. Finally, a mea- 
surement of the ancilla qubit produces |1) with probability 
C'^k{T). The precision of k{T) thus obtained goes as 1/M 
with the number of measurements if we use amplitude es- 
timation. The previously proposed approach for estimating 
reaction rates OlTIl evaluates the rate constant by computing a 
flux-flux correlation function based on a polynomial-size sam- 
ple of the wavefunction in the position basis. In contrast, our 
approach carries out the integration explicitly on the quan- 
tum computer using amplitude amplification, which provides 
a quadratic improvement over algorithms that rely on classical 
sampling and post-processing. 

The thermal state p(0) can be prepared efficiently on a 
quantum computer. We begin by preparing a superposi- 
tion of the reactant state labels in terms of C, and E, i.e., 
CEc,£;r(£^,T) IC-B), with C and V as defined above. 
Here, |C, £■) contain only the state labels, and not position- 
space wavefunctions. If we assume that the thermally acces- 
sible states can be enumerated by a polynomial number of 
qubits and that the energy can be specified to a certain pre- 
cision AE, we see that the states |C, £-) require only polyno- 
mially many qubits to store. The superposition itself can be 
prepared efficiently since T{E, T) is an exponential function 
of the energy and is therefore efficiently preparable [15, 31]. 

The next step is to generate, by doing state 
preparation in superposition, the state |<I>o) — 
cj2<:,E^iE,T)\C,E)\cj)o{C.E)), in which each term 
is the tensor product of \C,E) and the corresponding real- 
space reactant eigenstate |0o {C-. E)). The states |0o (d E)) 
must have definite energy E. Hence, one can represent an ini- 
tial state as a direct product of discrete reactant internal states 
(specified by ^ with energy E{C)) and a wavepacket with 
kinetic energy Ek — E — E{C) |39]. The discrete part can 
be prepared using the state-preparation approach described 
above. The incoming wavepacket can be approximated with 
a Gaussian with a kinetic energy expectation value of Ek- 
This approximation can be improved by increasing the width 
of the Gaussian, which can be doubled by the addition of a 
single qubit to the position register. This would halve the 
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momentum uncertainty of the wavepacket. With sufficiently 
many qubits, the error incurred in this approximation could 
be made smaller than errors from other sources, such as 
grid discretization. Once we have prepared \^o), we will no 
longer use the register containing the states \C,E). If we trace 
this register out, we can see that the remaining state register 
is a mixed state with density operator p{0), as desired. 

Finally, we show how to obtain state-to-state transition 
probabilities. Most chemical reactions can be regarded as 
scattering processes, and it is therefore desirable to know the 
scattering S'-matrix. In particular, it is these state-to-state tran- 
sition amplitudes which are accessible to experiment. Hereto- 
fore we have considered the joint wavefunction of all the 
molecular species. To compute state-to-state probabilities, 
however, we must ensure that each reactant molecule is in a 
well-defined molecular state. For example, to probe the state- 
to-state dynamics of the H + H2 reaction, we would need to 
prepare a particular state of the hydrogen atom plus a state of 
the hydrogen molecule, and not simply a state of the overall 
H3 aggregate. Given these states, prepared in the center-of- 
mass coordinate systems of each molecule, one must perform 
coordinate transformations to put them on a common grid. 
For each molecule, the Cartesian coordinates of the particles 
are linear, invertible functions of their center-of-mass coordi- 
nates. Since the coordinate transformations can be computed 
by an efficient, reversible, classical circuit, they can also be 
efficiently computed quantum mechanically. 

We concentrate here on obtaining only the vibrational state- 
to-state distributions. Using the techniques of state prepara- 
tion above, we prepare each reactant molecule in a vibrational 
eigenstate, so that along each of its normal mode coordinates 
the molecule is in an eigenstate of the corresponding potential. 
After all the molecules are thus prepared, their wavefunctions 
are transformed to a common Cartesian system. This initial 
state is evolved as usual until the molecules separate into iso- 
lated products. 

At large inter-molecular separation, the center-of-mass mo- 
tion and the normal mode coordinates again become separa- 
ble. Therefore, an orthogonal transformation can be applied 
to each product molecular fragment so that its Cartesian coor- 
dinates can be transformed into normal modes. The quantum 
phase estimation algorithm can then be employed to extract 
the populations and eigenenergies of the product vibrational 
states. 

For an isolated product molecule, we can expand the final 
state in terms of the normal modes: \'^') = J^v' <^v' l^'v')' 
where 1^'^/) is the position representation of the eigenstate 
corresponding to product occupation number vector v'. The 

state-to-state transition probabilities are then Pv'^v = lo^v'l^' 
and as mentined above, they can be determined using the 
phase estimation algorithm of Abrams and Lloyd 1 33] for each 
degree of freedom. We can obtain good measurement statis- 
tics with only a polynomial number of measurements because 
at typical temperatures, the products of chemical reactions 
will have appreciable population in only a small number of 
vibrational eigenstates. 



CONCLUSION 

The advantages of the methods presented here are not lim- 
ited to chemical reaction dynamics, but can be applied to 
many areas of physics. This is true in particular because the 
complexity of the algorithm is proportional to the complex- 
ity of calculating the potential and because the laws of nature 
are usually captured by simple, few-body interactions. For 
example, by using a quantum computer to study atoms in the 
presence of a strong, time-dependent electric field, one could 
simulate such effects as multielectron ionization or attosec- 
ond pulse generation 12811 . Quantum computers also 

offer the promise of predicting real-time properties of super- 
fluids Ii40l klj] , and of providing tests for effective potentials 
for water phases 14211 . 

We close by reiterating the need for a careful reexamination 
of the suitability of traditional quantum approximations for 
use on quantum computers. Previously we have shown that 
a quantum implementation of full configuration interaction 
scales better than coupled cluster approaches (in particular 
CCSD(T)), and in this work we show that simulating the com- 
plete nuclear and electronic time evolution is more efficient 
on quantum computers than using the Born-Oppenheimer ap- 
proximation, a central tool of theoretical chemistry. We can 
imagine the development of a wide variety of new tech- 
niques and approaches tailored for natural quantum simula- 
tors, which, themselves relying on the rules of quantum me- 
chanics, give us a deeper understanding of physical phenom- 
ena. 
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(project W9 11 NF-07- 1-0304 and the QuaCCR program) and 
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SUPPLEMENTARY INFORMATION: QUANTUM 
RESOURCE ESTIMATION 

There are two kinds of quantum register in the simulation 
algorithm: the state registers, used to store the wavefunction, 
and the ancilla registers, used to store the potential energy and 
the intermediate calculation results. If we assume a simulation 
of a d-dimensional system in which each Cartesian coordinate 
is divided into a uniform grid of A'' = 2" points, the represen- 
tation of the wavefunction requires a total of nd qubits in d 
registers. As for the ancilla registers, their total number will 
depend on the complexity of evaluating the potential and the 
kinetic energy. At least one register is always required, to be 
used as the target of addition for the purpose of phase kick- 
back. The ancilla registers will require m qubits each, where 
m is chosen in such a way that the registers can store the value 
of V{x) with desired accuracy, in the form of a binary integer 
between and 2"*. 

The time required for the simulation is the number of el- 
ementary (one- and two-qubit) gates required to perform the 
algorithm. Except in trivial cases, the evaluation of the po- 
tential energy will be much more complicated than that of the 
kinetic energy T, which is a simple quadratic form. We there- 
fore approximate the total gate count as being equal to the 
complexity of evaluating the potential: even for the simple 
Coulomb potential, the error thus introduced to the resource 
count is substantially less than one percent. 

Coulomb potential. The simulation of chemical dynamics 
depends on computing the Coulomb potential, and here we 
provide a detailed count of the resources required for evaluat- 
ing it on a quantum computer. We begin by developing some 
necessary quantum arithmetic. 

For addition, we adopt Draper's quantum addition algo- 
rithm [1], which is based on the quantum Fourier transform 
(QFT), and requires only controlled rotations. While 
it is not asymptotically optimal (i.e., it scales as 0{m?) and 
not 0{m) as does the schoolbook addition algorithm), it both 
has a small prefactor that makes it attractive for the addition 
of small numbers, and it is easily adapted for multiplication. 
Subtraction requires the same number of rotations, except that 
they are performed in the opposite direction. 

We perform multiplication using the schoolbook method. 
The first multiplicand is repeatedly bit shifted and added to 
the accumulator if the corresponding bit of the second multi- 
plicand is 1. Since each number has m bits, we need to make 
a total of m such controlled additions (C-ADD). The prod- 
uct will have 2m bits, but we will only keep the m most sig- 
nificant ones, essentially performing floating-point arithmetic. 
For the C-ADD, we first apply a QFT to the accumulator, as 
in Draper's algorithm (we will also apply an inverse QFT at 
the end, and these two require v? steps in all). Each C-ADD 
requires CC-rotations, which each can be implemented 
using two CNOTs and three C-rotations. Hence, each C-ADD 
requires operations, giving a total of -|- m^) gates 
for a multiplication. However, since half of the CC-rotations 
are to the insignificant bits of the accumulator and are sub- 
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sequently discarded, we only need to perform (|m'^ + m^) 
gates for a multiplication. 

To compute the Coulomb potential, the distance r between 
two particles, = (x2 - XiY + (2/2 - Vi)^ + (^2 - ziY, 
must be known. Evaluating requires 3 subtractions and 
3 squarings (the two additions are performed automatically 
since the squarings are really additions that can use the same 
accumulator). For squaring, we multiply the number by itself 
using the multiplication circuit, giving the total requirement 
of (^m^ + ^iTi^) gates for computing r^. The same com- 
putation would be used in momentum space for computing 
or for simulating a harmonic oscillator potential. 

The evaluation of the Coulomb potential is complicated by 
the need for a square root. Since evaluating is just as dif- 
ficult as evaluating 1/VS, we can find 1/r from in one 
computation using the Newton-Raphson method, with the it- 
eration Xn+i = ^ (3 — • x^). The number of iterations 
will depend on the desired final accuracy, but numerical ex- 
periments show that for many ranges of S, four iterations suf- 
fice to compute l/x/S* to within less than 0.03% over the en- 
tire range. Each iteration requires one subtraction and three 
multiplications (one of them bit-shifted due to the factor of 
^). So the requirement for l/\/S is (15m^ + 18m^) gates, 
which together with calculating the distance gives the to- 
tal requirement for the Coulomb potential as [-^m^ + ^rri^) 
gates for each pair of particles. 

Potentials fitted from first-principles calculations. When 
using the Born-Oppenheimer approximation, one uses a po- 
tential y(x) which is a function of only the nuclear coordi- 
nates. It is the total energy of the molecule assuming that the 
electrons are in their ground state given the potential induced 
by the nuclei at coordinates x. In general, this ground state 
energy is difficult to compute on a classical computer. Thus, 
interpolation schemes may be used to approximate V{x.). 
Here we analyze the computational resources needed for such 
schemes. 



We represent the potential as a d-dimensional interpolating 
polynomial: 

K 

ki ,k2 . ..kfi—0 
K K K 

fci=0 k2=0 kd=Q 

which can be evaluated using Horner's method starting with 
the innermost sum. That is, one has to evaluate one order- if 
polynomial in xi, K such polynomials in X2, and so on until 
K"^^^ polynomials in Xd- That is, the total number of polyno- 
mials that need to be evaluated is J2i=i ~ order 
that the results of these calculations be available as constants 
for higher-level polynomials, all the intermediate polynomial 
evaluations have to be saved in temporary memory along the 
way. The number of required registers is ■ 
Each polynomial that is evaluated has the form 

K 

P{x) = ^PkX^ = po+x{pi hx{pk-2 + X {pk-1 +Xpk))) 

k=0 

where pk is some constant. Therefore, each polynomial eval- 
uation requires precisely K additions and K multiplications. 
As we have seen above, an addition requires |to^ opera- 
tions, and a multiplication (|m'^ + m^), meaning that in total 
each polynomial evaluation requires K [jiri^ + |m^) gates. 
Therefore, the total number of gates required to calculate V 
is (fm^ + fm^) « K''- (|m^ + fm^). Furthermore, 

the total qubit requirement is ritot = nd + m ^_^ . 

1. T G. Draper, [quant-ph/000803 3 1 (2000) . 



